In den Geowissenschaften ist die Fernerkundung die einzige Messtechnik, die eine vollständige meßtechnische Abdeckung auch globaler Skalen ermöglicht. Zur erfolgreichen Nutzung gehören die souverände Anwendung von existierenden und notwendigerweise auch die Anpassung dieser bzw.die Entwicklung eigener Methoden.

Einleitung

Ein wichtiger Anwendungsfall in der Geo- oder Umweltinformatik ist die Erfassung von Veränderungen mittels Satelliten-, Flugzeug- und Drohnenbildern (Change Detectin Analysis). Oft wird dies in Verbindung mit bio- und geophysikalischen oder auch vom Menschen verursachten Prozessen genutzt, um ein tieferes Verständnis und die Möglichkeit zur Entwicklung von Vorhersagemodellen zu erhalten. Dabei sind Bildanalysemethoden von hervorgehobener Bedeutung, um Informationen zu extrahieren, die es erlauben, die zugrunde liegenden Prozesse zu identifizieren. Ein ständig wachsender Bereich ist die Integration von immensen Datenmengen sog. Big Data in die Analysen.

Die nahezu exponetiell wachsende und bereits gegenwärtig kaum noch zu bewältigende Flut an verfügbaren Bild- und Fernerkundungsdaten muss einfach und automatisiert zugänglich sein um sowohl für den wissenschaftlichen Erkenntnisgewinn als auch für gesellschaftliche Zukunftsaufgaben nutzbar zu sein. Wir beginnen mit einer anwendungs orientierten Aufgabe in dieser Übung.

Informationen aus Bilddaten

Unverarbeitete Satellitenbilder sind nicht zwingend informativ. Unser Auge kann zwar ein Echtfarbenbild relativ schlüssig und intuitiv interpretieren, aber eine zuverlässige und reproduzierbare wissenschaftliche Interpretation erfordert andere Ansätze. Zudem können Bild analysemethoden zusätzliche,unsichtbare und spezifischere Informationen ableiten. Wir haben bereits einfache Indizes berechnet. Auch die Ableitung des NDVI oder der Oberflächenalbedo ist eine physikalisch begründete Umwandlung von Bildsignalen in eine Meßgröße.

Um nützliche oder zielführende Informationen, z. B. über die Bodenbedeckung in einem Gebiet, zu erhalten, müssen wir die Daten daher fragestellungszentriert analysieren. Der wohl bekannteste und gängigste Ansatz ist die überwachte Klassifizierung von Bilddaten in Kategorien die von Interesse sind.

Diese Übung führt Sie in die Klassifizierung von Satelliten- und Luftvermessungsdaten in R ein.

Wir werden uns mit den folgenden Themen beschäftigen:

Vorbereiten der Arbeitsumgebung und Laden der Daten
Digtalisierung von Trainingsbereichen
Unüberwachte Klassifizierung (kmeans clustering)
Modelltraining
Überwachte Klassifizierung (random forest, Maximum Likelihood)
Güteschätzung Modell

Klassifizierung von Fernerkundungsdaten

Bitte beachten Sie, dass alle Arten der Klassifizierung eine ind der Regel aufwendige Datenvorverarbeitung erfordern. Im Zentrum stehen dann Modellbildung und Qualitätsschätzung die als handwerkliche Grundlagen der Klassifizierung betrachtet werden können um schließlich in dr Datennachverarbeitung als wesentlichen Bestandteil die inhaltliche Interpretation der Ergebnisse durchzuführen. Wir werden versuchen die notendigen Minimalanforderungen in der Kurseinheiten zu erarbeiten.

Überwachte Klassifizierung

Bei der überwachten Klassifizierung von Landbedeckungen wird aus einer begrenzten Menge Trainingsdaten zur Landbedeckung ein Modell abgeleitet, das die jeweilige Landbedeckung für den gesamten Datensatz vorhersagt. Die Landbedeckungstypen werden also a priori definiert, und das Modell versucht, diese Typen auf der Grundlage der Ähnlichkeit zwischen den Eigenschaften der Trainingsdaten und dem Rest des Datensatzes vorherzusagen.

Ganz pragmatisch erfordern Klassifizierungsaufgaben im Allgemeinen die folgenden Schritte:

Zusammenstellung eines umfassenden Eingabedatensatzes, der eine oder mehrere Rasterebenen enthält.
Auswahl von Trainingsgebieten, d.h. Teilmengen von Eingabedatensätzen, für die der Fernerkundungsexperte den Landbedeckungstyp kennt. Das Wissen über die Landbedeckung kann z.B. aus eigenen oder fremden in situ Beobachtungen, Managementinformationen oder anderen Fernerkundungsprodukten (z.B. hochauflösenden Luftbildern) gewonnen werden.
Training eines Modells unter Verwendung der Trainingsflächen. Zu Validierungszwecken werden die Trainingsflächen häufig in eine oder mehrere Test- und Trainingsstichproben unterteilt, um die Leistung des Modellalgorithmus zu bewerten.
Anwendung des trainierten Modells auf den gesamten Datensatz, d. h. Vorhersage der Bodenbedeckungsart auf der Grundlage der Ähnlichkeit der Daten an jedem Ort mit den Klasseneigenschaften des Trainingsdatensatzes.

Die folgende Abbildung zeigt die Schritte einer überwachten Klassifikation im Detail. Die optionalen Segmentierungsoperationen sind obligatorisch für objektorientierte Klassifizierungen, die sich nicht nur auf die Geometrie der Objekte, sondern auch auf die Werte jeder einzelnen Rasterzelle im Eingabedatensatz stützen. Um einzelne Objekte wie Häuser oder Baumkronen abzugrenzen, verwenden Fernerkundungsexperten Segmentierungsalgorithmen, die die Homogenität der Pixelwerte innerhalb ihrer räumlichen Nachbarschaft berücksichtigen.

Change Detection Waldveränderung Nord-West-Harz

In diesem Tutorium werden die Sentinel-2-Bilder aus der vorherigen Übung verwendet.

Start - Einrichten des der Arbeitsumgebung

Sie können entweder die gespeicherten Daten aus der vorangegangenen Einheit verwenden oder einen neuen Raumausschnitt bestimmen, herunterladen und bearbeiten. Im Prinzip wird jedoch zuerst die Arbeitsumgebung geladen.

#  ---- 0 Projekt Setup ----
# Achtung Pakete müssen evtl. manuell installiert werden
library(envimaR)

#--- Schalter für den Download der sentinel daten
get_sen = FALSE

#--- schalter ob digitalisiert werden muss falls er auf FALSE gesetzt ist werden die
# (zuvor erstellten und gesciherten Daten ) im else Teil der Verzeigung eingelesen
digitize = FALSE

## setzen des aktuellen Projektverzeichnisses (erstellt mit envimaR) als root_folder
#root_folder = find_rstudio_root_file()
root_folder = "~/edu/geoinfo/"

# Einlesen des zuvor erstellten Setup-Skripts
source(file.path(root_folder, "src/functions/000_setup.R"))

nclasses=2

Bitte ergänzen Sie (bei auftreteneden Fehlermeldungen) etwaig fehlende oder defekte Pakete im obigen Setup-Skripts.

Auf der Grundlage der verfügbaren Sentinel Daten sollten nun als Erstes geeignete Datensätze für eine Oberflächenklassifikation identifiziert werden. Hierzu kann der vollständige Datensatz auch vom Kursdatenserver heruntergeladen werden (Bitte beachten Sie dass sie im VPN bzw. UniNetz angemeldet sein müssen). Entpacken Sie diese Daten in das Wurzelverzeichnis des Kursprojekts. Das heisst der dort bereits vorhandene Ordner data wird ersetzt/ergänzt.

Schritt 1 Übersicht verschaffen

Bei näherer Betrachtung der RGB Bilder (RGB432B) zeigt sich, das vier Datensätze aufgrund der Bildqualität und geringen Wolkenbedeckung geeignet zu sein scheinen. Es handelt sich um den 19.06. bzw. 24.7. 2019 und den 33.06 bzw. 30.7. 2020. Aufgrund des früheren Vegetationszeitpunktes wurden die Maibilder gewählt, da hier evtl. Aufwuchs auf den Rodungsflächen weniger stark sichtbar sein könnte.

Zunächst müssen diese Daten in einem “Rasterstapel” also einem Mehrkanalbild verfügbar sein. Bereits bei der ersten Beschäftigung mit dem sen2r Paket hatten wir weitere Produkte als sinnvoll erachtet und berechnet. Es sind folgende Indizes:

Normalized Difference Vegetation Index NDVI
Modified Soil Adjusted Vegetation Index MSAVI2
Simple Ratio 1600/820 Moisture Stress Index MSI
Modified Soil Adjusted Vegetation Index SAVI
Enhanced Vegetation Index EVI

Diese Daten müssen für den gewünschten Zeitslot und Raumausschnitt besorgt und vorbereitet werden. Aus Gründen der Vereinfachung gibt es für den Kurs eine entsprechende sen2R-Projektdatei:

{
  "preprocess": [true],
  "s2_levels": ["l1c", "l2a"],
  "sel_sensor": ["s2a", "s2b"],
  "online": [true],
  "server": ["scihub"],
  "order_lta": [true],
  "downloader": ["builtin"],
  "overwrite_safe": [false],
  "rm_safe": ["yes"],
  "max_cloud_safe": [5],
  "step_atmcorr": ["auto"],
  "sen2cor_use_dem": [true],
  "sen2cor_gipp": [null],
  "timewindow": ["2019-06-01", "2020-08-31"],
  "timeperiod": ["seasonal"],
  "extent": ["{\n  \"type\": \"FeatureCollection\",\n  \"features\": [\n    {\n      \"type\": \"Feature\",\n      \"properties\": {\n        \"X_leaflet_id\": 2080,\n        \"layerId\": \"2080\",\n        \"edit_id\": \"2080\"\n      },\n      \"geometry\": {\n        \"type\": \"Polygon\",\n        \"coordinates\": [\n          [\n            [10.198975, 51.848575],\n            [10.198975, 51.947721],\n            [10.413895, 51.947721],\n            [10.413895, 51.848575],\n            [10.198975, 51.848575]\n          ]\n        ]\n      }\n    }\n  ]\n}"],
  "s2tiles_selected": ["32UNC"],
  "s2orbits_selected": [null],
  "list_prods": ["BOA", "SCL", "CLD"],
  "list_indices": ["ARVI", "EVI", "MSAVI2", "MSI", "NDVI", "SAVI"],
  "list_rgb": ["RGB432B", "RGB843B"],
  "rgb_ranges": [
    [null],
    [
      [0, 7500],
      [0, 2500],
      [0, 2500]
    ]
  ],
  "index_source": ["BOA"],
  "mask_type": ["nodata"],
  "max_mask": [5],
  "mask_smooth": [0],
  "mask_buffer": [0],
  "clip_on_extent": [true],
  "extent_as_mask": [true],
  "extent_name": ["harz"],
  "reference_path": [null],
  "res": [null],
  "res_s2": ["10m"],
  "unit": ["Meter"],
  "proj": [null],
  "resampling": ["near"],
  "resampling_scl": ["near"],
  "outformat": ["GTiff"],
  "rgb_outformat": ["GTiff"],
  "index_datatype": ["Int16"],
  "compression": ["DEFLATE"],
  "rgb_compression": ["DEFLATE"],
  "overwrite": [true],
  "path_l1c": ["~/edu/geoinfo/data/data_lev0"],
  "path_l2a": ["~/edu/geoinfo/data/data_lev0"],
  "path_tiles": [null],
  "path_merged": [null],
  "path_out": ["~/edu/geoinfo/data/data_lev1"],
  "path_rgb": ["~/edu/geoinfo/data/data_lev1"],
  "path_indices": [""],
  "path_subdirs": [true],
  "thumbnails": [true],
  "log": ["~/edu/geoinfo//doc/nw-harz.log"],
  "parallel": [true],
  "processing_order": ["by_date"],
  "pkg_version": ["1.5.0.9000"]
}

#--- Download der Daten
# gui = TRUE ruft die GUI zur Kontrolle auf
if (get_sen){
   out_paths_3 <- sen2r(
    gui = T,
    param_list = "~/edu/geoinfo/data/harz.json",
    tmpdir = envrmt$path_tmp,
  )
}
#--- Einlesen der Daten aus den Verzeichnissen
# RGB stack der beiden Jahre

pred_stack_2019 = raster::stack(list.files(file.path(envrmt$path_data_lev1,"RGB843B"),pattern = "20190619",full.names = TRUE))
pred_stack_2020 = raster::stack(list.files(file.path(envrmt$path_data_lev1,"RGB843B"),pattern = "20200623",full.names = TRUE))

# Stack-Loop über die Daten
for (pat in c("RGB432B","EVI","MSAVI2","NDVI","SAVI")){
  pred_stack_2019 = raster::stack(pred_stack_2019,raster::stack(list.files(file.path(envrmt$path_data_lev1,pat),pattern = "20190619",full.names = TRUE)))
  pred_stack_2020 = raster::stack(pred_stack_2020,raster::stack(list.files(file.path(envrmt$path_data_lev1,pat),pattern = "20200623",full.names = TRUE)))
}

# Zuweisen von leserlichen Namen auf die Datenebenen
names(pred_stack_2019) = c("nir","red","green","red","green","blue","EVI","MSAVI2","NDVI","SAVI")
names(pred_stack_2020) = c("nir","red","green","red","green","blue","EVI","MSAVI2","NDVI","SAVI")
saveRDS(pred_stack_2019,paste0(envrmt$path_data,"pred_stack_2019.rds"))
saveRDS(pred_stack_2020,paste0(envrmt$path_data,"pred_stack_2020.rds"))
pred_stack_2019 = readRDS(paste0(envrmt$path_data,"pred_stack_2019.rds"))
pred_stack_2020 = readRDS(paste0(envrmt$path_data,"pred_stack_2020.rds"))
names(pred_stack_2019) = c("nir","red_1","green_1","red_2","green_2","blue","EVI","MSAVI2","NDVI","SAVI")
names(pred_stack_2020) = c("nir","red_1","green_1","red_2","green_2","blue","EVI","MSAVI2","NDVI","SAVI")

# visuelle Überprüfung der stacks
plot(pred_stack_2019)

plot(pred_stack_2020)

# Einlesen der Corine Daten
# Für den Download ist ein Konto notwendig https://land.copernicus.eu/pan-european/corine-land-cover
# Daher die Daten manuell herunterladen und in das Verzeichnis kopieren und entpacken
# Dann auskommentierten Tail ausführen
# corine_eu = raster(file.path(envrmt$path_data_lev0,"u2018_clc2018_v2020_20u1_raster100m/DATA/U2018_CLC2018_V2020_20u1.tif"))
# tmp = projectRaster(pred_stack_2019[[1]],crs = crs(corine_eu))
# corine_crop = raster::crop(corine_eu,tmp)
# corine_utm = projectRaster(corine_crop,crs = crs(pred_stack_2019))
# corine = resample(corine_utm,pred_stack_2019[[1]])
# raster::writeRaster(corine,file.path(envrmt$path_data_lev0,"/corine.tif"),overwrite=TRUE)

# ansonsten den Beipieldatensatz laden
corine = raster::raster(file.path(envrmt$path_data_lev0,"corine.tif"))

# Erstellen einer Wald-Maske
# Agro-forestry areas code=22, Broad-leaved forest code=23,
# Coniferous forest code=24, Mixed forest code=25
mask = reclassify(corine,c(-100,22,0,22,26,1,26,500,0))

mapview(mask)

## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; 
## the supplied Raster* has 1693870 
##  ... decreasing Raster* resolution to 5e+05 pixels
##  to view full resolution set 'maxpixels =  1693870 '

K-Means Cluster-Klassifikation

Die wohl bekannteste unüberwachte Klassifizierungstechnik ist das K-means-Clustering, das auch als der “einfachster Algorithmus des maschinellen Lernens” bezeichnet wird. Er wird häufig angewendet um einen ersten Überblick zu erhaltne ob die Rasterdaten im Merkmalsraum ausreichend trennbar sind.

In unserem Beispiel (auf 2 Klassen angewendet und mit unsuperClass Funktion aus dem RStoolbox Paket ausgeführt) sieht das wie folgt aus:

## k-means über RStoolbox
# Modell
prediction_kmeans_2019 = unsuperClass(pred_stack_2019, nClasses = 2,norm = TRUE, algorithm = "MacQueen")
# Klassifikation
mapview(prediction_kmeans_2019$map, col = c('darkgreen', 'burlywood', 'green'))

## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; 
## the supplied Raster* has 1693870 
##  ... decreasing Raster* resolution to 5e+05 pixels
##  to view full resolution set 'maxpixels =  1693870 '

prediction_kmeans_2020 = unsuperClass(pred_stack_2020, nClasses = 2,norm = TRUE, algorithm = "MacQueen")
mapview(prediction_kmeans_2020$map, col = c('darkgreen', 'burlywood', 'green'))

## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; 
## the supplied Raster* has 1693870 
##  ... decreasing Raster* resolution to 5e+05 pixels
##  to view full resolution set 'maxpixels =  1693870 '

Schritt 2 - Trainingsdaten erstellen

Trainingsdaten werden vielfach durch abdigitalisieren manuell erfasst. Dies kann recht komfortabel in RStudio durchgeführt werden, falls nur schnell und effektiv einige wenige Trainingsflächen zu digitalisieren sind.

Für größere Aufgaben ist es sinnvoll, auf den hohen Komfort z.b. der QGIS-Arbeitsumgebung zurückzugreifen. Für diese Übung verwenden wir mapedit, ein kleines, aber feines Paket, das die Digitalisieren und Editieren von Vektordaten im Rstudio- oder externen -browser ermöglicht. In Kombination mit mapview können auch beliebige Farbkomposita als Grundlage für die Digitalisierung verwendet werden.

Trainingsdaten digitalisieren

Wir nehmen an, dass wir zwei Arten von Landbedeckung klassifizieren wollen: clearcut (Abholzungen) und other (Anderes). Mit mapedit `muss jede Klasse einzeln digitalisiert werden. Sobald die Trainingsgebiete als Vektordaten verfügbar sind können die Merkmale des jeweiligen Rasterstacks entsprechend der digitalisierten Klassen in eine Tabelle extrahiert und auf etwaige Fehlwerte bereiningt werden.

Falls dieser Teil bereits absolviert wurde kann die logische Variable (am Anfang des Scripts definiert) digitize auf FALSE gesetzt werden und es wird dann der else Teil der Verzweigung durchlaufen - also nur noch die existierenden Daten eingelesen.

Verwendung von Farbkomposita für bessere Trainingsergebnisse

m1 = tm_shape(pred_stack_2019) + tm_rgb(r=1, g=2, b=3) +
  tm_layout(legend.outside.position = "right",
            legend.outside = T,
            panel.label.height=0.6,
            panel.label.size=0.6,
            panel.labels = c("r=1, g=2, b=3")) +
  tm_grid()

m2 = tm_shape(pred_stack_2019) + tm_rgb(r=4, g=5, b=6) +
  tm_layout(legend.outside.position = "right",
            legend.outside = T,
            panel.label.height=0.6,
            panel.label.size=0.6,
            panel.labels = c("r=4, g=5, b=6")) +
  tm_grid()
tmap::tmap_arrange(m1,m2)

## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)

## Warning: Raster values found that are outside the range [0, 255]

## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)

## Warning: Raster values found that are outside the range [0, 255]

## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)

## Warning: Raster values found that are outside the range [0, 255]

## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)

## Warning: Raster values found that are outside the range [0, 255]

Verwenden Sie die Ebenensteuerung, um die Ebenen umzuschalten. Bei Echtfarbkompositen werden die sichtbaren Spektralkanäle Rot (B04), Grün (B03) und Blau (B02) den entsprechenden roten, grünen bzw. blauen Farbkanälen zugeordnet, wodurch ein quasi natürliches “farbiges” Bild der Oberfläche entsteht, wie es ein Mensch sehen würde, der auf dem Satelliten sitzt. Falschfarbenbilder werden häufig mit den Spektralkanälen für das nahe Infrarot, Rot und Grün erzeugt. Sie eignen sich hervorragend für die Einschätzung der Vegetation, da Pflanzen nahes Infrarot und grünes Licht reflektieren, während sie rotes Licht absorbieren (Red Edge Effect). Ein dichterer Pflanzenbewuchs ist dunkler rot. Städte und offener Boden sind grau oder hellbraun, und Wasser erscheint blau oder schwarz.

#---- Digitalisierung der Trainingsdaten ----

if (digitize) {
  # Für die überwachte Klassifikation benötigen wir Trainingsgebiete. Sie können Sie wie nachfolgend digitalisieren oder alternativ z.B. QGis verwenden
  
  #--- 2019
  # Kahlschlag
  # Es wird das Falschfarbenkomosit in originaler Auflösung genutzt (maxpixels =  1693870)
  # Bitte beachten Sie dass es (1) deutlich länger lädt und (2) Vegetation in Rot dargestellt wird.
  # Die Kahlschäge sind jetzt grün
  train_area <- mapview::viewRGB(pred_stack_2019, r = 1, g = 2, b = 3, maxpixels =  1693870) %>% mapedit::editMap()
  # Hinzufügen der Attribute class (text) und id (integer)
  clearcut <- train_area$finished$geometry %>% st_sf() %>% mutate(class = "clearcut", id = 1)
  
  # other: hier gilt es möglichst verteilt übers Bild möglichst alle nicht zu Kahlschlag  gehörenden Flächen zu erfassen.
  train_area <- mapview::viewRGB(pred_stack_2019, r = 1, g = 2, b = 3) %>% mapedit::editMap()
  other <- train_area$finished$geometry %>% st_sf() %>% mutate(class = "other", id = 2)
  
  # rbind  kopiert die beiden obigen Vektorobjekte in eine Datei
  train_areas_2019 <- rbind(clearcut, other)
  
  # Umprojizieren auf die Raster Datei
  train_areas_2019 = sf::st_transform(train_areas_2019,crs = sf::st_crs(pred_stack_2019))
  mapview(train_areas_2019, zcol="class")
  
  # Extraktion der Trainingsdaten für die digitalisierten Flächen
  tDF_2019 = exactextractr::exact_extract(pred_stack_2019, train_areas_2019,  force_df = TRUE,
                                          include_cell = TRUE,include_xy = TRUE,full_colnames = TRUE,include_cols = "class")
  #  auch hier wieder zusamenkopieren in eine Datei
  tDF_2019 = dplyr::bind_rows(tDF_2019)
  
  # Löschen von etwaigen Zeilen die NA (no data) Werte enthalten
  tDF_2019 = tDF_2019[complete.cases(tDF_2019) ,]
  tDF_2019 = tDF_2019[ ,rowSums(is.na(tDF_2019)) == 0]
  
  # check der extrahierten Daten
  summary(tDF_2019)
  mapview(train_areas_2019)+pred_stack_2019[[1]]
  
  # Abspeichern als R-internes Datenformat
  # ist im Repo hinterlegt und kann desahlb (zeile drunter) eingeladen werden
  saveRDS(tDF_2019, paste0(envrmt$path_data,"train_areas_2019.rds"))
  
  
  # # ---- Das gleiche muss für 2020 wiederholt werden zum digitalisieren und extrahieren bitte ent-kommentieren ----
  
  # Kahlschlag
  train_area <- mapview::viewRGB(pred_stack_2020, r = 1, g = 2, b = 3,maxpixels =  1693870) %>% mapedit::editMap()
  clearcut <- train_area$finished$geometry %>% st_sf() %>% mutate(class = "clearcut", id = 1)
  train_area <- mapview::viewRGB(pred_stack_2020, r = 1, g = 2, b = 3) %>% mapedit::editMap()
  other <- train_area$finished$geometry %>% st_sf() %>% mutate(class = "other", id = 2)
  train_areas_2020 <- rbind(clearcut, other)
  train_areas_2020 = sf::st_transform(train_areas_2020,crs = sf::st_crs(pred_stack_2020))
  mapview(train_areas_2020, zcol="class")
  tDF_2020 = exactextractr::exact_extract(pred_stack_2020, train_areas_2020,  force_df = TRUE,
                                          include_cell = TRUE,include_xy = TRUE,full_colnames = TRUE,include_cols = "class")
  tDF_2020 = dplyr::bind_rows(tDF_2020)
  tDF_2020 = tDF_2020[  rowSums(is.na(tDF_2020)) == 0,]
  saveRDS(tDF_2020, paste0(envrmt$path_data,"train_areas_2020.rds"))
  
} else {
  train_areas_2019 = readRDS(paste0(envrmt$path_data,"train_areas_2019.rds"))
  train_areas_2020 = readRDS(paste0(envrmt$path_data,"train_areas_2020.rds"))
  
}

Als Resultat ligen nun zwei tabellarische Trainingsdatensätze für 2019 und 2020 vor. Beide Datensätze beinhalten alle Rasterinformationen die von den Polygonen für die Klassen “clearcut” und “other” abgedeckt werden.

head(train_areas_2019)

##      class nir red_1 green_1 red_2 green_2 blue  EVI MSAVI2 NDVI SAVI      X
## 1 clearcut  74    64      58    64      58   39 2969   2648 5524 2978 586485
## 2 clearcut  67    55      49    55      49   32 2778   2498 5684 2854 586495
## 3 clearcut  75    96      73    96      73   49 2191   2018 3975 2300 586475
## 4 clearcut  75   109      78   109      78   57 1962   1792 3464 2049 586485
## 5 clearcut  74    97      73    97      73   53 2161   1954 3876 2234 586495
## 6 clearcut  62    72      56    72      56   38 2100   1897 4422 2221 586505
##         Y   cell coverage_fraction
## 1 5752735 514566        0.27757052
## 2 5752735 514567        0.08858948
## 3 5752725 516064        0.32311472
## 4 5752725 516065        0.99449599
## 5 5752725 516066        0.98001540
## 6 5752725 516067        0.68315113

head(train_areas_2020)

##      class nir red_1 green_1 red_2 green_2 blue  EVI MSAVI2 NDVI SAVI      X
## 1 clearcut  86    95      72    95      72   53 2793   2542 4606 2820 586645
## 2 clearcut  83   102      74   102      74   52 2475   2297 4190 2569 586655
## 3 clearcut  57    83      54    83      54   39 1576   1453 3473 1732 586665
## 4 clearcut  85    99      72    99      72   50 2603   2427 4390 2701 586645
## 5 clearcut  85   100      71   100      71   51 2615   2435 4389 2708 586655
## 6 clearcut  75    99      68    99      68   48 2154   2007 3915 2285 586665
##         Y   cell coverage_fraction
## 1 5753385 417147       0.249708712
## 2 5753385 417148       0.474895865
## 3 5753385 417149       0.006561131
## 4 5753375 418646       0.750996768
## 5 5753375 418647       1.000000000
## 6 5753375 418648       0.717114925

Schritt 3 - Modelltraining und Prüfen der Modellgüte

Klassifikatoren (z.B der Maximum-Likelihood Klassikator) oder Algorithmen des maschinellen Lernens (wie z. B. Random Forest) ermitteln auf Basis der Trainingsdaten Beschreibungsmodelle die z.b. statistische Signaturen, Klassifikationsbäume oder andere Funktionen darstellen. Im Rahmen der Güte der Trainingsdaten sind solche Modelle geeignet und repräsentativ um Vorhersagen für Räume zu treffen falls die Prädiktoren aus dem Modell flächendeckend vorhanden sind.

Wir wollen nun die räumlichen Merkmale Kahlschlag/kein Wald exemplarisch mit einer Maximum Likelihood Klassifikation und mit Random Forest vorhersagen und Standardmethoden der Zufallsvalidierung und Modellgüteeinschätzung anwenden.

Ziel ist es Kahlschläge von allen anderen Pixeln zu trennen und die Unterschide von 2019/2020 zu quantifizieren.

Random forest

Random Forests können sowohl für Regressions- als auch für Klassifizierungsaufgaben verwendet werden, wobei letztere besonders in der Umwelt-Fernerkundung relevant sind. Wie bei jeder maschinellen Lernmethode lernt auch das Random-Forest-Modell, Muster und Strukturen in den Daten selbst zu erkennen.

Abbildung: Vereinfachte Darstellung der Klassifizierung von Daten durch Random Forest während des Trainings. Venkata Jagannath [CC BY-SA 4.0] via wikipedia.org

Ein Random-Forest-Algorithmus lernt über die Daten, indem er viele Entscheidungsbäume erstellt - daher auch der Name “Forest.” Für Klassifizierungsaufgaben nimmt der Algorithmus eine passende Instanz aus eines Baumes aus dem Trainingsdatensatz und weist dem Pixel dies Klasse zu. Dies wird mit allen verfügbaren Entschidungsbäuumen so gemacht und letzlich wird nach dem Winner takes it all Ansatz die Instanz der Klasse zugeordnet, die das Ergebnis der meisten Bäume ist.

Da der Random-Forest-Algorithmus Trainingsdaten benötigt, handelt es sich um eine überwachte Lernmethode. Das bedeutet, dass wir als Nutzer dem Algorithmus Daten zur Verfügung stellen müssen, die im Wissen über die vorherzusagenden Klassen vermitteln. Diese Daten werden dann in Trainings- und Testdaten aufgeteilt.

## ---- Überwachte  mit Klassifikation Random Forest ----
  train_areas_2019 = readRDS(paste0(envrmt$path_data,"train_areas_2019.rds"))
  train_areas_2020 = readRDS(paste0(envrmt$path_data,"train_areas_2020.rds"))
## Hier wird der Random Forest über das Utility Paket caret aufgerufen
# Setzen eines "seed" ermöglicht reproduzierbaren Zufall
set.seed(123)

# Zufälliges Ziehen von 25% der Daten (Training/Test)
idx_2019 = createDataPartition(train_areas_2019$class,list = FALSE,p = 0.25)
trainDat_2019 = train_areas_2019[idx_2019,]
testDat_2019 = train_areas_2019[-idx_2019,]
idx_2020 = createDataPartition(train_areas_2020$class,list = FALSE,p = 0.25)
trainDat_2020 = train_areas_2020[idx_2020,]
testDat_2020 = train_areas_2020[-idx_2020,]

#  Response-Variable (=Spalte "class") muss den Datentyp "factor" haben
trainDat_2019$class <- as.factor(trainDat_2019$class)
testDat_2019$class <- as.factor(testDat_2019$class)
trainDat_2020$class <- as.factor(trainDat_2020$class)
testDat_2020$class <- as.factor(testDat_2020$class)


# Einstellungen Modelltraining: cross-validation, 10 Wiederholungen
ctrlh = trainControl(method = "cv",
                     number = 10,
                     savePredictions = TRUE)

#--- random forest model training
cv_model_2019 = train(trainDat_2019[,2:11], # in den Spalten 2 bis 20 stehen die Trainingsdaten (Prediktoren genannt)
                      trainDat_2019[,1],         # in der Spalte 1 steht die zu klassizierende Variable (Response genannt)
                      method = "rf",             # Methode hier rf für random forest
                      metric = "Kappa",          # Qualitäts/Performanzmaß KAppa
                      trControl = ctrlh,         # obig erzeugte Trainingssteuerung soll eingelsen werden
                      importance = TRUE)         # Die Bedeung der Variablen wird mit abgespeichert

cv_model_2019

## Random Forest 
## 
## 5996 samples
##   10 predictor
##    2 classes: 'clearcut', 'other' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 5397, 5396, 5397, 5396, 5396, 5396, ... 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##    2    0.9879911  0.8578546
##    6    0.9884903  0.8658856
##   10    0.9879903  0.8600513
## 
## Kappa was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 6.

Berechnung einer Wahrheitsmatrix zur Überprüfung der Modellgüte

Die Testdaten werden nun für die unabhängige Qualitätsprüfung des Modells verwendet. Eine Wahrheits- oder Konfusionsmatrix zeigt an, wie präzise das Modell die korrekten Klassen vorhersagt. Die Hauptdiagonale der Matrix zeigt die Fälle an, in denen das Modell korrekt ist. In diesem Falle mit nur 2 Klassen gilt der Sonderfall der Beurteilung eines binären Klassifikators. Speziell für die hier verwendete Funktion finden sich ausführliche Erläuertungen in der caret Hilfe.

Die wesentlichen Ausagen zur Modellgüte sind : * ‘Positive’ Class = clearcut: wird mit der Sensitivität (true positive rate) erfasst, die die Wahrscheinlichkeit angibt, mit der ein positives Objekt korrekt als positiv klassifiziert wird. * ‘Negative Class’ = other: wird mit der Spezifität (true negative rate) erfasst und gibt die Wahrscheinlichkeit an, mit der ein negatives Objekt korrekt als negativ klassifiziert wird * Positive and negative predictive values geben für clearcut bzw. other die reale Perfomance an. Sie sind um die reale Häufikeitsverteilung korrigiert und ein Schätzmass für die Präzision bzw. Pereformance des Modells bezüglich der jeweiligen Klassen.

Trotz der hohen Werte sehen wir das die Klasse clearcut hier deutlich abfällt also durchaus noch Verbesserungsbedarf für eine guten Klassifikation vorliegt.

Insgesamt kann das Modell jedoch als gut bezeichnet werden.

# ---- Berechnung der Konfusionsmatrix  ----
cm_rf <- confusionMatrix(data = predict(cv_model_2019, newdata = testDat_2019), testDat_2019$class)
cm_rf

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction clearcut other
##   clearcut      699   103
##   other         120 17063
##                                           
##                Accuracy : 0.9876          
##                  95% CI : (0.9859, 0.9892)
##     No Information Rate : 0.9545          
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0.8559          
##                                           
##  Mcnemar's Test P-Value : 0.284           
##                                           
##             Sensitivity : 0.85348         
##             Specificity : 0.99400         
##          Pos Pred Value : 0.87157         
##          Neg Pred Value : 0.99302         
##              Prevalence : 0.04554         
##          Detection Rate : 0.03887         
##    Detection Prevalence : 0.04459         
##       Balanced Accuracy : 0.92374         
##                                           
##        'Positive' Class : clearcut        
##

Vorhersage auf die Ausgangsdaten

Nun sind wir soweit das überprüfte Modell auf unseren Datensatz anzuwenden. Üblicherweise wird das in der Fernerkundung Klassifikation genannt.

# Klassifikation (auch Vorhersage genannt)
prediction_rf_2019  = raster::predict(pred_stack_2019 ,cv_model_2019 )

#--- 2020
# Einstellungen wie 2019
cv_model_2020 = train(trainDat_2020[,2:11],
                      trainDat_2020[,1],
                      method = "rf",
                      metric = "Kappa",
                      trControl = ctrlh,
                      importance = TRUE)

cv_model_2020

## Random Forest 
## 
## 7845 samples
##   10 predictor
##    2 classes: 'clearcut', 'other' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 7061, 7060, 7060, 7060, 7061, 7060, ... 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##    2    0.9847046  0.9218005
##    6    0.9841944  0.9192613
##   10    0.9838122  0.9174582
## 
## Kappa was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 2.

prediction_rf_2020  = raster::predict(pred_stack_2020 ,cv_model_2020)

Maximum Likelihood Klassifikation

Bei der Maximum-Likelihood-Klassifizierung wird davon ausgegangen, dass die Verteilung der Daten für jede Klasse und in jedem Kanal normal verteilt sind. Unter dieser Vorraussetzung wird die Wahrscheinlichkeit berechnet, dass ein bestimmtes Pixel zu einer bestimmten Klasse gehört. Da auch die Wahrscheinlichkeiten als Schwellenwert angegeben werden können werden ohne diese Einschränkung alle Pixel ungeachtet wie unwahrscheinlich zugeordnet. Jedes Pixel wird der Klasse zugeordnet, die die höchste Wahrscheinlichkeit aufweist (d. h. die maximale Wahrscheinlichkeit).

# ---- Maximum Likelihood Classification ----
# superClass() Funktion aus dem Paket RSToolbox umwandeln der Tabelle in das
# geforderte SpatialdataPoint Objekt Identifikation der Koordinaten
sp::coordinates(trainDat_2019) = ~X+Y
sp::coordinates(trainDat_2020) = ~X+Y
# Zuweisen der korrekten Projektion (hier aus dem zugehörigen Raster-Stack)
crs(trainDat_2019) = crs(pred_stack_2019)
crs(trainDat_2020) = crs(pred_stack_2020)

# superClass method "mlc" berechnet die Statisik und klassifiziert in einem aufruf
prediction_mlc_2019       <- superClass(pred_stack_2019, trainData = trainDat_2019[,1:11], responseCol = "class",
                                        model = "mlc", tuneLength = 1, trainPartition = 0.5,verbose = FALSE)

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

prediction_mlc_2020       <- superClass(pred_stack_2020, trainData = trainDat_2020[,1:11], responseCol = "class",
                                        model = "mlc", tuneLength = 1, trainPartition = 0.5,verbose = FALSE)

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

# Ergebnisse der rf Klassifikation
#prediction_mlc_2019$model
#prediction_mlc_2020$model


## ---- Visualisierung mit mapview ----
mapview::viewRGB(mask*pred_stack_2020, r = 4, g =5, b = 6,maxpixels =  1693870)+
  mapview(mask*prediction_rf_2019 , alpha.regions = 0.5, maxpixels =  1693870,
          col.regions = mapviewPalette("mapviewRasterColors"),at = seq(0, nclasses, 1), legend = TRUE) +
  mapview(mask*prediction_rf_2020, alpha.regions = 0.5, maxpixels =  1693870,
          col.regions = mapviewPalette("mapviewRasterColors"),at = seq(0, nclasses, 1), legend = FALSE) +
  mapview(mask*prediction_mlc_2019$map,alpha.regions = 0.5, maxpixels =  1693870,
          col.regions = mapviewPalette("mapviewRasterColors"),at = seq(0, nclasses, 1), legend = FALSE) +
  mapview(mask*prediction_mlc_2020$map,alpha.regions = 0.5, maxpixels =  1693870,
          col.regions = mapviewPalette("mapviewRasterColors"),at = seq(0, nclasses, 1), legend = FALSE)

## Warning in raster::projectRaster(x, raster::projectExtent(x, crs =
## sp::CRS(epsg3857)), : input and ouput crs are the same

Ressourcen

Betrachten sie sie nachfolgenden Ressourcen als Beispiele dafür, wie aus solchen Quellen (die alle mehr oder weniger technisch ähnlich sind) ein bestimmtes Set von Werkzeugen zur Bearbeitung einer Fragestellung “entsteht”. Dann nach viel Recherche und kritischer Prüfung wird ein aktueller Stand der Technik innerhalb der wissenschaftlcihen Gemeinschaft erreicht der sich zudem natürlich stetig weiterentwickelt.

Schauen sie einmal in die nachstehende Auswahl von Blogs und Hilfestellungen:

Digitalisierung von Trainingsdaten by the Forest Inventory and Remote Sensing department at the University of Goettingen (Germany)
Digitalisierung Turorial in the QGIS 3.16 documentation
Robert J. Hijmans rspatial - supervised classification
Ivan Lizarazo RPubs Tutorial
Sydney Goldstein blog
João Gonçalves supervised classification
Valentin Stefan pixel-based supervised classification

Die beiden ersten Beiträge sind rein technische Tutorials. Die weiteren Beiträge mischen technische Anleitungen mit konzeptionellen oder konkreten fachlichen Fragestellungen. Dennoch sind sie hier nicht als fachwissenschaftliche Referenz gedacht. Sie dienen vielmehr dazu zu zeigen wie technisches und konzeptionelles Verständnis schrittweise erarbeiet werden kann und unterstützen, durch “nachkochen” und anwenden, eigenständiger die Lösung von Fragestellungen anzugehen.

Ich möchte ausdrücklich Valentin Stefan den Autor des Blogbeitrags pixel-based supervised classification zitieren:

“[…] Betrachten Sie diesen Inhalt als einen Blogbeitrag und nichts weiter. Er erhebt nicht den Anspruch, eine erschöpfende Übung oder ein Ersatz für Ihr kritisches Denken zu sein. […]”